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Abstract 

Quantitative real time RT-PCR has been described as the most sensitive method for the detection of low abundance mRNA. 
To date, no reference genes have been screened in the half-smooth tongue sole (Cynoglossus semilaevis). The aim of this 
study was to select the most stable genes for quantitative real-time RT-PCR. Eight housekeeping genes (18S, TUBA, B2M, 
ACTB, EF1A, GAPDH, RPL17 and UBCE) were tested at different developmental stages, in different tissues, and following 
exposure to the drug SB-431542. Using geNorm, BestKeeper and NormFinder software, GAPDH/B2M, GAPDH/18S and 
UBCE/GAPDH were identified as the most suitable genes from samples taken of different developmental stages while 18S/ 
RPL17 were consistently ranked as the best reference genes for different tissue types. Furthermore, TUBA/B2M, TUBA/UBCE 
and B2M/TUBA were found to be the most suitable genes in samples treated with the drug, SB-431542 by geNorm, 
BestKeeper and NormFinder respectively. Across both different developmental stages and tissue types, the combination of 
18S and GAPDH was the most stable reference gene analyzed by Ref-Finder. To test and verify the screened reference 
genes, the expression profiles of LEFTY-normalized to the combination of GAPDH/18S and ACTB were presented. These 
results will be useful for future gene-expression studies in the half-smooth tongue sole. 
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Introduction 

Quantitative real-time PCR (qRT-PCR) plays an extremely 
important role in studies on gene expression and determination of 
complex molecular pathways in various biological systems [1]. It is 
a robust method with obvious advantages, including sensitivity, 
large dynamic range, and the potential for high throughput and 
accurate quantification [2] . The read-out can be either an absolute 
number of copies, or a relative amount normalized to an internal 
control gene [2]. The accuracy of this technique, however, might 
be affected at multiple stages throughout the experimental process 
and by several factors. It has been proposed that housekeeping 
genes might be used as reference genes for normalizing the results 
of qRT-PCR because of their consistent expression, both with 
respect to different developmental stages and exposure to different 
environmental factors. In reality, however, no single reference 
gene can be used for longitudinal studies for different tissues, 
ontogenetic stages and experimental conditions. The transcript 
levels of genes deemed as housekeeping traditionally may also vary 
considerably [3]. Thus, the optimal genes can only be confirmed 
with an analysis of models based on particular cases. An accepted 



method used to reduce errors to a minimum is to normalize the 
target gene to a set of reference genes identified through an 
analysis of parameter models [1,4]. GeNorm, BestKeeper and 
NormFinder are Visual Basic Applications that were developed for 
running in Microsoft Excel and they are widely used in various 
conditional experiments [4]. Pairwise variation is calculated by 
geNorm to determine the most stable reference genes, whereas 
BestKeeper is used to evaluate the coefficient of variation. 
NormFinder takes the analysis of intra- and inter-group variations 
for normalization into consideration. 

The half-smooth tongue sole (Cynoglossus semilaevis) is an 
important farmed marine fish, which shows behavioral and 
morphological left-right asymmetry, with one eye migrating to 
the opposite side (metamorphosis) during the larval stages [5]. A 
high larval mortality due to starvation can occur during 
metamorphosis when body malformations cause the mouth to fail 
to close. Whereas the mechanism regulating this external 
asymmetry remains largely unknown, evidence suggests that the 
Nodal signaling pathway is involved in the eye laterality of a 
flatfish similar to the sole, the Japanese flounder [6-9]. More 
studies are needed to analyze the expression profiles of genes 
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related to the Nodal signaling pathway. Currently, litde work has 
been conducted on the selection of reference genes in the half- 
smooth tongue sole, although some reference genes have been 
identified in the Atlantic halibut and Japanese flounder [10-12]. 
While it is useful to consider the reference genes identified from 
the halibut and flounder, inter-species variation can occur and so it 
is more accurate to identify the particular genes for the species 
under investigation, in this case the half-smooth tongue sole [13]. 

The aim of this work was to screen optimal sets of reference 
genes for 18 different developmental stages, 8 tissue types and 
following treatment with SB-431542, an inhibitor of the Nodal 
signaling pathway. Three statistical algorithms, geNorm, Best- 
Keeper and NormFinder, were used to analyze the qRT-PCR 
data. In addition, the ontogenesis expression profiles of LEFTY 
were investigated to confirm the efficacy of the selected reference 
genes as this is reported to be part of the Nodal-pathway that plays 
a role in controlling eye laterality in the flounder [7], as well as 
being shown to be expressed on the left side of the dorsal 
diencephalon and internal organs to form fixed laterality [6,8]. 
The results obtained in this study provide an essential tool for 
future gene expression studies in the half-smooth tongue sole. 

Materials and Methods 

Ethics Statement 

The half-smooth tongue sole used in this study were collected 
from local aquatic farms under a permit from the local 
government of Yantai, Shandong, China. All handling of the 
half-smooth tongue sole was conducted in accordance with the 
guidelines and regulations established by the Ocean University of 
China and the local government of YanTai. 

Fish Rearing and Sample Collection 

Fertilized eggs were obtained by artificial fertilization and then 
maintained at 22°C. The developmental stages were determined 
by observation via microscopy. A total of 1 8 developmental stages 
were selected: 1 -cell (0.5 hours post fertilization; hpf), 4-cell 
(1.5 hpf), 8-cell (1.8 hpf), 16-cell (2.3 hpf), 32-cell (2.8 hpf), 1 It- 
cell (6 hpf) and sphere (11 hpf) stages, as well as 30% epiboly 
(13.5 hpf), and the 2-somite (20 hpf), 15-somite (3.5 hpf), 21- 
somite (27.5 hpf), 27-somite (30 hpf) and metamorphosis stages D, 
E, F, G, H, and I, as described by Minami (1982) [14]. In brief, the 
metamorphosis stages can be described as follows: Stage D (the 
stage prior to the start of eye migration, 15 DAH; day after 
hatching); Stage E (the eye begins to migrate, 1 7 DAH); Stage F 
(the migrating eye is visible from the ocular side, 18 DAH); Stage 
G (the upper edge of the migrating eye is beyond the dorsal 
margin, 19 DAH); Stage H (the upper edge of the migrating eye is 
beyond the dorsal mid-line, 20 DAH); Stage I (the entire 
migrating eye is past the dorsal mid-line, 22 DAH) (Minami, 
1982). The larvae were raised at 22°C. 

These 1 8 stages were arranged into 3 groups each containing 6 
stages, called developmental stage I (from the 1 -cell to lk-cell 
stage), II (sphere stage to 27-somite stage) and III (metamorphosis 
stages D-I). For developmental stages I and II, 30 embryos were 
used per stage, while for developmental stage III, three larvae were 
used. In chemical treated experiments, embryos were treated from 
the 16-cell stage to the 21 -somite stage with SB-431542 at 100 uM 
and 200 uM (or with double-distilled water or DMSO as controls) 
with 30 embryos being treated per experimental group. All 
experiments were performed in biological triplicate. 

Regarding the different tissues, a total of six adult fish (i.e., three 
female and three male fish) were involved in this investigation. 
Eggs and sperm came from three female and three male fish, 



respectively, while the other tissues (muscle, kidney, gill, spleen, 
liver, intestines, egg and sperm) were selected from 3 different 
adult fish randomly, and again experiments were performed in 
triplicate, and then stored in liquid nitrogen immediately until use. 
All procedures complied with the Institutional Animal Care and 
Use Committee of Ocean University of China and were 
specifically approved for this study. 

Selection of Reference Genes 

A total of eight candidate genes were selected for the gene 
expression stability analysis as follows: 18S (Ribosomal RNA), 
TUBA (oc-Tubulin), ACTB (P-actin), B2M (fS-2-Microglobulin), 
EF1A (Elongation factor-l-a), GAPDH (Glyceraldehyde-3-phos- 
phate), RPL17 (Ribosomal protein LI 7) and UBCE (Ubiquitin- 
conjugating enzyme). Moreover, LEFTY (Fw: 5'- 
CCC AAGGTGTTGTC C GTCTG-3 ' ; Rev: 5'-ATT- 
CAAGGCTGATCCCC ATTC-3') was used to verify the proper 
reference genes. All the genes were identified from an unpublished 
cDNA library of Cynogbssus semilaevis. The primers were designed 
with Primer Premier 5.0 software (PREMIER Biosoft Interna- 
tional, Palo Alto, CA). The full gene name, function, primers and 
PCR parameters are shown in Table I. The length of PCR 
products ranged from 70 to 170 bps, and the mean melting 
temperature of the forward and reverse primers ranged between 
61.0°C and 62.7°C. 

RNA Extraction, cDNA Synthesis and RT-PCR 

Total RNA was extracted from samples using Trizol reagent 
according to the manufacturer's instructions. The concentration 
and purity of the total RNA extracted were determined using a 
Nano Photometer. The absorbance ratio of each sample at OD 
260/280 was between 1.9 and 2.0, and the OD 260/230 was 
approximately 2.0. The integrity of total RNA was assessed via 1 % 
ethidium bromide agarose gel electrophoresis. A total of 1 ug 
RNA from each sample was reverse transcribed according to the 
instructions of the PrimeScript™ RT reagent kit with gDNA 
Eraser (Takara, Japan), and the final volume was 20 ul. Standard 
RT-PCR was carried out using Premix Ex Taq to characterize all 
the primer pairs according to the length and specificity of the 
products. 

Quantitative Real-time PCR (qRT-PCR) 

Each 20 ul reaction volume contained 10 ul 2X SYBR Premix 
Ex Taq™, 0.4 ul ROX Reference Dye II, 0.4 ul of each primer 
(10 uM), 2 ul 10X diluted cDNA template (1000X diluted cDNA 
for 18S rRNA) and 6.8 ul nuclease-free water. The PCR 
conditions were as follows: 95°C for 30s followed by 40 cycles of 
95°C for 5s and then 60°C for 40s. PCR was performed using an 
Applied Biosystems 7500 Real-Time PCR system. Three technical 
replicates were performed for each sample. 

Data Analysis 

To comprehensively and systematically assess the expression 
variation of the candidate reference genes, all the samples were 
divided into three broad categories: developmental stages, 
different tissues, and chemical treatment. The developmental 
stages were subdivided into three groups: Stage I, the 1-cell stage 
(i.e., a newly fertilized egg) to the sphere stage; Stage II, 30% 
epiboly to the 27-somite stage; and Stage III, throughout 
metamorphosis. 

The time sequence fluorescence values (Rn) were collected using 
the qRT-PCR machine. All raw data, which based on all cDNAs 
samples of tissue, developmental stage, and chemical-treated, were 
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Table 1. Reference genes for gene expression normalization in Cynoglossus semilaevis. 





Gene 


Molecular function 


Primer (5 -3 ) 


Product 


E 


R 2 


Average 

T m (°C) 


symbol 






size (bp) 








18S 


Ribosomal subunit 


Fw: GGTAACGGGGAATCAGGGT 


70 


0.98 


0.989 


61.2 


Rv: TGCCTTCCTTGGATGTGGT 


TUBA 


Cytoskeletal protein, 


Fw: ATCCTCCTGCTGTTGCCTTTG 


107 


0.99 


0.990 


62.0 




microtubule activity 


Rv: TGACCTGATGTCGTATCGCCTC 










ACTB 


Cytoskeletal structural 


Fw: GAAATCGCCGCACTTGTTGT 


111 


1.00 


0.987 


61.6 




protein 


Rv: GGGTCAGGATACCTCTCTTGCTCT 










B2M 


Subunit of the MHC class 


Fw: TGTTCGTCGTTCTGCCGTGT 


112 


0.97 


0.990 


62.2 




1 molecule 


Rv: TCAGGGTGTTGGGCTTGTTGT 










EF1A 


Protein synthesis 


Fw: AGGCTGGTATCTCCAAGAACGG 


122 


0.98 


0.987 


62.5 


Rv: TCTGGCTGTAAGGGGGCTCA 


GAPDH 


Glycolytic enzyme 


Fw: GAAGGGCATTCTGGGATACACT 


119 


1.03 


0.993 


61.0 




phosphorylating activity 


Rv: TCAAAGATGGAGGAGCGGC 










RPL17 


Structural constituent of 


Fw: AGTGCGTCCCGTTCCGTC 


114 


1.00 


0.991 


61.5 




ribosomal subunit 


Rv: TCAGCGTTGCTCTCTGCGTT 










UBCE 


Ligase activity, protein 


Fw: GCAATGGTAGCATCTGCCTCG 


170 


1.00 


0.989 


62.7 




degradation 


Rv:GCTTCTTCGTTATACTTCTGTCGGTC 











Note: E means reaction efficiencies and R means Pearson's coefficients of determination. Average Tm represents the mean melting temperature of the forward (Fw) 
and reverse (Rv) primers. Each sample contains 3 biological replicates. 
doi:1 0.1 371 /journal.pone.0091 71 5.t001 



imported into the PCR Miner software to calculate efficiency (E), 
cycle threshold (Ct), and Pearson's coefficients of determination 
(R 2 ) of each reference gene. The expression stability level was 
determined with geNorm and NormFinder, which use relative 
expression values as input data by transformed the Ct values to 
linear scale expression quantities using the 2" delta " Ct method, and 
BestKeeper based on. the raw Ct value directly using the 
comparative ACT method. 

The overall final ranking of reference genes across all 
developmental stages and tissue types was calculated using Ref- 
Finder, which integrates the current major computational 
programs (geNorm, Normfinder, BestKeeper, and the compara- 
tive delt Ct method) to compare and rank the stability of candidate 
reference genes with the geometric mean of individual gene 
appropriate weight [15]. 

Results 

Expression Profiles of Candidate Reference Genes 

Eight candidate reference genes, including 18S, TUBA, B2M, 
ACTB, EF1A, GAPDH, RPL17 and UBCE, were selected for the 
evaluation of expression stability with n = 3 experiments contain- 
ing 30 samples per experiment for developmental stages I and II, 
and 3 samples per experiment for developmental stage III, as 
described in the Materials and Methods. The amplification 
specificity of the primers used in the qRT-PCR was evaluated 
using agarose gel electrophoresis followed by ethidium bromide 
(EB) staining. The specificity of the primers producing single 
amplicons was confirmed by the absence of primer dimers and 
false priming. Specific PCR products of identical length were 
identified via electrophoresis (Figure 1) and sequencing results. 
RT-PCR data analysis showed that the amplification efficiency (E) 
of the eight reference genes ranged from 0.97 for B2M to 1.03 for 
GAPDH, the mean coefficient variation (R 2 ) ranged from 0.9877 



for B2M to 0.9938 for GAPDH, and the mean melting 
temperature (Tm) of the forward and reverse primers ranged 
from 61.0°C for GAPDH to 62.7°C for UBCE (Table 1). 

Throughout all developmental stages (Figure 2A), the genes 
with the lowest global variability were GAPDH and B2M. 
GAPDH, ACTB and B2M were the top three genes with the 
lowest variation in developmental stage I (Figure 2B); 18S, 
GAPDH and B2M showed littie variation in stage II (Figure 2C 
and Figure A in File SI) and stage III (Figure 2D and Figure B in 
File SI). However, 18S, RPL17 and EF1A showed smaller 
variations in different tissues (Figure 2E and Figure C in File SI), 
and TUBA, B2M and UBCE were the least variable genes in the 
chemically treated group (Figure 2F and Figure D in File SI). 

To assess the Ct value of eight candidate reference genes, we 
divided the raw expression levels of the genes into three parts 
according the average Ct values: (1) Genes with a high transcript 
abundance exhibited an average Ct value below 20; (2) those with 
median transcript abundance had an average Ct value of 20-25; 
while (3) those with low transcript abundance had an average Ct 
value of 25-35. EF1A showed high transcript abundance in all 
groups except in stage I (all Ct values of eight reference genes were 
great than 20 in stage I); 18S showed median transcript 
abundance, while B2M, GAPDH showed low transcript abun- 
dance in all sample pools (Table SI). The expression level 
variations of the eight candidate genes indicated that optimized 
reference genes could be identified in each group for different 
research purposes. 

geNorm Analysis 

The geNorm analysis was used to select the best reference genes. 
Two parameters were considered to quantify reference gene 
stability: M value (average expression stability) and Vn/n+1 
(pairwise variation). Originally, MS1.5 and Vn/n+l<0.15 are 
regarded as acceptable levels of expression variability [4]. 
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lOObp 




Figure 1. Specificity of qRT-PCR and amplicon length. Amplified fragments shown on 1% agarose gel electrophoresis with ethidium bromide 
staining. An identical length of specific PCR products was demonstrated. M: Marker DL2000; Lanes 1-8:18S, TUBA, B2M, ACTB, EF1A, GAPDH, RPL17 
and UBCE, respectively. 
doi:1 0.1 371 /journal.pone.0091 71 5.g001 



Recently, an M value below the threshold of 0.5 is generally 
regarded as being typical for a stable reference gene in a relatively 
homogeneous sample panel. For more heterogeneous panels, the 
mean M value may increase to 1 [16]. Depending on the 
developmental stage, tissue type and chemical treatment, the 
housekeeping gene expression stability showed some variation. 
Throughout all developmental stages, the most stable genes were 
GAPDH and B2M. The M value obtained for these two genes was 
0.57, which is below the default limit of Msl. 18S was also 
identified as stable, with an M value of 0.66 (Figure 3A). The two 
most stable genes in developmental stage I were 18S and UBCE 
(Figure 3B), GAPDH and B2M in stage II (Figure 3C); GAPDH 
and ACTB in stage III (metamorphic stages) (Figure 3D), 18S and 
RPL17 for different tissues (Figure 3E) and B2M and TUBA for 
the chemical treatment group (Figure 3F). To determine the 
number of required reference genes for each group, the Vn/ n+ 1 
was evaluated using geNorm. All Vn/n+1 values, including V2/3, 
were below 0.15 in each pool (Figure 3G), which indicated that 
only two genes are needed as an internal control. 

BestKeeper Analysis 

Gene expression variation was calculated for all eight candidate 
reference genes based on Ct-values and displayed as the standard 
deviation (SD) and coefficient of variance (CV) by BestKeeper. 
The lowest values of SD (i.e., usually <1) indicated the highest 
stability [17]. BestKeeper highlighted GAPDH as the reference 
gene with the least overall variation from the list of eight candidate 
genes, with an SD of 0.29, CV value = 0.92 throughout all 
developmental stages (Figure 4A). The variation in expression of 
the other candidate reference genes was greater except the 18S 
and B2M (SD = 0.37, and SD = 0.38 respectively). The order of 
the candidate reference genes from most stable (lowest SD) to least 
stable (highest SD) was GAPDH >18S>B2M>UBCE>EF1A> 
AC TB > TUB A> RPL 1 7 . Following for stage I, the ranking was 
GAPDH>B2M>ACTB >18S>UBCE>EF1A>RPL17> TU- 
BA (Figure 4B); 18S>GAPDH>B2M>TUBA>EF1A>UBCE> 
RPL17> ACTB for stage II (Figure 4C); and 18S>GAPDH> 
UBCE>B2M>ACTB>EF1 A>TUBA>RPL1 7 for metamorphic 
stage III (Figure 4D). For different tissue types, the decreasing 
order of stability was 18S>RPL17> EF 1 A> GAPDH > ACTB > 
TUBA>UBCE>B2M (Figure 4E). 18S (SD = 0.45, CV= 1.86) 
showed constant expression stability (Table S2); The order of 
variation changed to UBCE>TUBA>RPL1 7> B2M>EF1A > 
18S>ACTB>GAPDH in the chemically treated group 
(Figure 4F). The most stable gene was UBCE (SD = 0.26, 



CV= 1.02) (Table S2), and the top five genes had similar stability 
levels, with SD values ranging from 0.26 to 0.4. 

NormFinder Analysis 

NormFinder was applied to validate the most stable expression 
gene [18]. As shown in Table 2, a high M value represents a high 
expression variance. Thus, UBCE, GAPDH and 18S occupied the 
top 3 positions for the most stable expression genes in the group 
throughout developmental stages I, II and III. The most suitable 
combination of genes was UBCE and GAPDH (combined stability 
value = 0.632). The rest were ranked as Mows: TUBA, B2M, 
EF1A, RPL 17 and ACTB. Based on developmental stages, the 
ranking was 18S>GAPDH>UBCE>EF1A>ACTB>RPL17> 
B2M>TUBA for stage I. The most suitable combination of genes 
was 18S and GAPDH (combined stability value = 0.154); 
GAPDH>EF1A >18S>B2M>UBCE>TUBA>RPL17> 
ACTB for stage II. The most suitable combination of genes was 
GAPDH and EF1A (combined stability value = 0.203); and 
GAPDH >18S>B2M>ACTB>UBCE>TUBA>RPL17> 
EF1A for stage III. The most suitable combination of genes was 
GAPDH and 18S (combined stability value = 0.038). For different 
tissue types, the sequence of stability was RPL17>18S>EF1A> 
ACTB>GAPDH>TUBA>UBCE>B2M. The most stable com- 
bination of genes was RPL 17 and 18S (combined stability 
value = 0.264). The ranking order of stability changed to B2M> 
TUBA>EF1A>UBCE >18S>RPL17> GAPDH>ACTB in the 
chemically-treated group. The most stable combination of genes 
were B2M and TUBA with the combined stability value = 0.039. 

According to the comprehensive ranking of the Ref-Finder 
performed, the combination of 18S and GAPDH with geomean 
1.19 and 1.73 were confirmed as the most stable genes other than 
EF1A (2.62), RPL 17 (4.76), UBCE (5.23), B2M (5.66), TUBA 
(5.73) and ACTB (7.74) across all developmental stages and tissue 
types. 

Expression Profile of LEFTY Gene of Ontogenesis and 
Drug Treatment 

To test and verify the results of the reference gene screening, 
which were analyzed by geNorm, Bestkeeper and NormFinder, 
expression profiling of LEFTY was carried out via qRT-PCR. 
Four parameters, GAPDH, 18S, a normalization factor based on 
the geometric average of both Ct values (GAPDH and 18S) and 
ACTB were used to normalize the LEFTY qRT-PCR data. They 
showed a similar tendency during the continuous development 
process: a low and stable level from fertilization to the 32-cell 
stage, after which the level increased more or less consistently from 
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Figure 2. The relative level of mRNA expression of the eight candidate reference genes. The data are presented as mean cycle threshold 
(Ct) values and shown as box-whisker plots. The boxes represent the inter quartile range of the mean Ct values, while the bars represent the 
maximum and minimum Ct values. The dotted line represents the mid-hinge and the scattered points represent the abnormal data. (A) 
Developmental stages I - III; (B) developmental stage I; (C) developmental stage II; (D) metamorphosis stage III; (E) different tissues; (F) chemical 
treatment. 

doi:1 0.1 371 /journal.pone.0091 71 5.g002 



high stage, to peak at the 21 -somite stage, followed by a reduction 
in expression from the 27-somite stage to approximately 1 day 
after hatching, after which low-to-negligible levels were observed 
during the metamorphosis stages. However, the expression level 
identified by the combination of GAPDH andl8S was more 
reliable when gene expression showed different trends between 
single reference gene group and the combination group. For 
example, LEFTY gene expression showed no significant decrease 
from 21 somites to 27 somites labeled by GAPDH, but it 



significantly dropped by the combination of GAPDH and 18S. 
The expression result revealed by the combination group was 
more consistent with the report in Japanese flounder. In addition, 
the expression level of LEFTY was lower compared with 
GAPDH/ 18S, and either GAPDH or 18S when ACTB was as 
reference gene through the entirely developmental stages. For 
instance, at 21 somites stage, no significant expression difference of 
LEFTY was shown among GAPDH, 18S and GAPDH/ 18S 
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Figure 3. Expression stability (M) analysis and the optimal reference gene number determination with geNorm. The stability values of 
the eight candidate genes at: (A) Developmental stages l-lll; (B) developmental stage I; (C) developmental stage II; (D) metamorphosis stage III; (E) 
different tissues; (F) chemical treatment; (G) The optimal number of reference genes determined with geNorm, whereby Vn was calculated from 
leastwise 2 genes. Vn/n+1 present the pairwise variation between 2 sequential Vn and Vn+1. 
doi:1 0.1 371 /journal.pone.0091 71 5.g003 



groups with high transcript abundance, but significantly lower 
expression was found when ACTB as reference gene (Figure 5). 

Discussion 

The qRT-PCR method is widely used in gene expression studies 
of both freshwater [19-21] and marine fish [22,23], when 
investigating a variety of conditions, such as infection [24,25] 
and stress [20,26-28], as well as in various cell lines [29-31] and 
different tissues [21-23,32]. As a crucial step, the qRT-PCR data 
should be accurately normalized by use of the appropriate 
reference genes. Otherwise, the use of non-validated reference 
genes can lead to erroneous conclusions that are biologically 
meaningless [33]. However, there is no universally applicable 
reference gene. The same reference genes in different species and 
different reference genes in the same species can have dramatically 
different expression levels under diverse experimental conditions. 
For example, GAPDH has been validated as the most stable 
reference gene in the mandarin fish [34], but has significant 
expression variation in the Atlantic halibut [33] and zebrafish [35]. 
On the other hand, ACTB and Tubb2 showed less expression 
variation, but 18S rRNA and GAPDH varied significandy, even 
these reference genes in the same species, Atlantic halibut [36]. 
Thus, to improve the sensitivity and reliability of qRT-PCR in 



practice, the validation of reference genes should be performed for 
each different species, tissue type and experiment condition to be 
used. 

The results obtained in our study were analyzed with geNorm, 
BestKeeper and NormFinder. The most stable reference genes 
generated from these three software applications showed a high 
level of similarity; however there were subde deviations in the rank 
order. For example, throughout all developmental stages, the top 3 
stable expression genes were: GAPDH, B2M and 18S for geNorm; 
GAPDH, 18S and B2M for Bestkeeper; and 18S, GAPDH and 
UBCE for NormFinder. This deviation is natural because of the 
different algorithms used. Both BestKeeper and geNorm use a 
pair-wise comparison approach, and are highly dependent on the 
assumption that none of the genes being analyzed are co-regulated 
[36]. On the other hand, in NormFinder, a separate analysis of the 
sample subgroups and estimation of both intra- and inter-group 
variation in expression levels are included into the calculation of a 
gene stability value [33]. 

In our investigation, GAPDH and B2M, were identified as the 
most suitable genes in samples containing all eighteen develop- 
mental stages, with an M-value of 0.57 (i.e., M < 1); while 18S was 
the next most stable gene with an M-value of 0.66 (i.e., M<1). In 
order to obtain a more accurate view about the stability of 
reference genes for the different developmental stages, we sub- 
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Figure 4. Expression stability analysis of the 8 candidate reference genes by BestKeeper. The stability values of the eight candidate genes 
at: (A) Developmental stages l-lll; (B) developmental stage I; (C) developmental stage II; (D) metamorphosis stage III; (E) different tissues; (F) chemical 
treatment. 

doi:1 0.1 371 /journal.pone.0091 71 5.g004 



divided the original eighteen developmental stages into three 
groups: termed I, II and III and we showed that the optimum 
reference genes in each group possessed M-values of less than 0.5 
(Figure 3B-D). To ensure that the normalization is as accurate as 
possible, reference genes with the closest Ct value to that of the 



target genes should be selected [20] . Thus, the expression level of 
each reference gene was also investigated (Table SI). 

According to the comprehensive ranking of the Ref-Finder 
performed, the combination of two reference genes, 18s and 
GAPDH, were shown to be the most stable genes across all 



Table 2. Expression stability analysis of the 8 candidate reference genes by NormFinder. 





Rank 


Stages l-lll 


Stage 1 


Stage II 


Stage III 


Tissues 


Treatment 


1 


UBCE 


185 


GAPDH 


GAPDH 


RPL17 


B2M 


M value 


0.778 


0.193 


0.196 


0.059 


0.240 


0.039 


2 


GAPDH 


GAPDH 


EF1A 


18S 


18S 


TUBA 


M value 


1.372 


0.278 


0.481 


0.104 


0.288 


0.039 


3 


18S 


UBCE 


18S 


B2M 


EF1A 


EF1A 


M value 


1.380 


0.467 


0.816 


0.164 


0.754 


0.123 


4 


TUBA 


EF1A 


B2M 


ACTB 


ACTB 


UBCE 


M value 


1.564 


0.542 


0.993 


0.178 


1.760 


0.139 


5 


B2M 


ACTB 


UBCE 


UBCE 


GAPDH 


18S 


M value 


1.827 


0.712 


1.008 


0.203 


1.973 


0.22 


6 


EF1A 


RPL17 


TUBA 


TUBA 


TUBA 


RPL17 


M value 


1.964 


0.792 


1.034 


0.255 


3.193 


0.269 


7 


RPL17 


B2M 


RPL17 


RPL17 


UBCE 


GAPDH 


M value 


2.077 


0.900 


1.250 


0.283 


3.355 


1.710 


8 


ACTB 


TUBA 


ACTB 


EF1A 


B2M 


ACTB 


M value 


2.376 


1.123 


1.958 


0.765 


4.331 


1.823 


the best 


UBCE/ 


18S/GAPDH 


GAPDH/ 


GAPDH/ 


RPL1 71 


B2M/ 


combination 


GAPDH 


0.154 


EF1A 


18S 


18S 


TUBA 


M value 


0.632 




0.203 


0.038 


0.264 


0.039 



Note: The stability of the eight candidate genes, indicated by the M value were calculated by NormFinder at all developmental stages (i.e., l-lll combined); as well as at 
developmental stage I; stage II; and stage III alone; in different tissues; and following chemical treatment. 
doi:1 0.1 371 /journal.pone.0091 71 5.T.002 
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Figure 5. Quantitative analysis of the LEFTY gene expression in embryos. The relative expression variance is shown as the ratio between 
the amount of LEFTY mRNA normalized to the corresponding GAPDH values, 18S value and the geometric average of both values at a particular 
developmental stage and in unfertilized eggs. Different developmental stages from 1 cell stage to metamorphosis stages. Columns correspond to: 1- 
cell, 4-cell, 8-cell, and 32-cell stages, high blastula stage, sphere stage, 30% epiboly, as well as the 2-somite, 15-somite, 21-somite, and 27-somite 
stage, during hatching, in 1-day old larva, and during metamorphosis stages D~l. 
doi:1 0.1 371 /journal.pone.0091 71 5.g005 



developmental stages and tissue types in the half-smooth tongue 
sole. This is consistent with the results reported in the Atlantic 
salmon [37] and Japanese flounder [12]. In addition, tissue- 
dependent variations have been observed in the expression of most 
housekeeping genes in flatfish. In our study, 18S and RPL17 were 
shown to be the most stable reference genes. ACTB has been used 
previously for half-smooth tongue sole organs [38-40], but it was 
proved to be an unstable reference gene in our investigation 
because of significant difference between ACTB and the 
combination of GAPDH and 18S, even either of 18S or GAPDH. 
ACTB was ranked at the worst reference gene in tissues group. In 
contrast, ACTB/UBCE were shown to be the best reference genes 
in different organs of the Japanese flounder [10], and EF1A and 
RPSD were identified for Atlantic halibut [36] and turbot [41] 
tissue. 

The expression stage of LEFTY was detected by in situ 
hybridization from shield-stage to the 27-somite stage in the 
Japanese flounder [42], We demonstrated that LEFTY appeared 
earlier (at high blastula, i.e., Ik-cell stage) in the half smooth 
tongue sole, and up to ldpf by qRT-PCR, which is a more 
sensitive assay than in situ hybridization. In zebrafish, there are two 
LEFTY genes (LEFTY 1 and LEFTY2). Most reports describe the 
expression of LEFTY 1 occurring from prior to gastrulation (from 
sphere stage in the blastula period) [43], to the 25- or 26-somite 
stage in the segmentation period [44-49] although there is no 
evidence showing the expression of either LEFTY 1 or LEFTY2 
after this stage. However, our evidence indicates that in the half 
smooth tongue sole, the expression of LEFTY appears earlier and 
disappears later (i.e., it is expressed over a wider range of 
developmental stages), when compared with either the zebrafish or 
the Japanese flounder. 



We believe that this is the first work to assess valid reference 
genes for qRT-PCR studies in the half-smooth tongue sole. This 
work is important for future developmental gene expression studies 
in this commercially important species as it provides valuable tools 
for investigating gene expression in both the embryonic and larval 
stages, as well as in different tissues and following chemical 
treatment in this flatfish species, which currently suffers from a 
high mortality risk during larval production. 

Supporting Information 

Figure SI The transcriptional levels (Ct value) of 8 
candidate reference genes (18S, TUBA, B2M, ACTB, 
EF1A, GAPDH, RPL17 and UBCE). Presented as mean Ct 
value (cycle threshold value) for each candidate gene and each 
sample of different developmental stages, tissues and chemical 
treatment. For developmental stages I and II, 30 embryos were 
used per stage, while for developmental stage III, three larvae were 
used. Tissue samples were collected from three different adult fish. 
Experiments were performed in triplicate. (A) Embryonic 
developmental stages (I and II), include the 1 -cell, 4-cell, 8-cell, 
16-cell, 32-cell, high-, and sphere-stages, 30% epiboly, and 2- 
somite, 15-somite, 21-somite, and 27-somite stages; (B) Metamor- 
phosis stages (III), includes metamorphosis stages D, E, F, G, H, I; 
(C) Different tissues include muscle, kidney, gill, spleen, liver, 
intestines, egg and sperm; (D) Embryos at the 16-cell stage were 
treated with SB431542 at 100 uM or 200 uM, or with DMSO or 
double distilled water as controls. 
(TIF) 

Table SI Groups divided according to abundance of 
expression of eight candidate genes. 
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Table S2 Descriptive statistics of 8 candidate reference 
genes expression based on the BestKeeper approach. 

(DOG) 
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